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Abstract 

Within the hquid drop model, and using a phenomenological Skyrme-Uke parameterization of 
the uniform nuclear matter equation of state (EoS), we first construct a set of 'baseline' crustal 
equations of state that span the experimentally constrained values of the symmetry energy at 
nuclear saturation density, as well as the theoretically constrained region of the pure neutron matter 
(PNM) EoS at low density. The resulting range of crustal compositional parameters relevant for 
macroscopic models of the crust is explored, as well as the crust-core transition densities and 
pressures and the density range of the 'pasta' phases of nuclei at the bottom of the crust. We 
explore the deviation of these parameters from the baseline models as a result of the remaining 
experimental uncertainties in the symmetric nuclear matter (SNM) EoS and the parameterization 
of the surface energy in the liquid drop model. The transition densities and pressures are found 
to be quite sensitive to the behavior of the surface tension at very low proton fractions; recent 
calculations of the energies of neutron drops suggest that this might be higher than previously 
thought, which our study suggests may result in a much reduced volume of pasta in the crust. 
Most crustal compositional parameters are shown to be insensitive to the SNM EoS and the surface 
energy relative to the uncertainties in the symmetry energy, with the notable exception of the size 
of the nuclei which show a similar dependence on all the model parameters. We establish a set 
of crustal equations of state which can be matched with any core EoS based upon the symmetry 
energy and its derivative, allowing a consistent exploration of the dependence on the symmetry 
energy of global crustal models and potentially observable phenomena such as pulsar glitches and 
oscillation modes that incorporate a coupling between crust and core. 

PACS numbers: 21.65.Cd,21.65.Ef,21.65.Mn,26.60.Gj,26.60.Kp,97.60.Jd 



I. INTRODUCTION 



The increasing wealth and sensitivity of observations associated with neutron star sys- 
tems is demanding ever more sophisticated theoretical models of their structure and dy- 
namics from micro- to macro-scopic scales. Many proposed theoretical models for observed 
phenomena, such as pulsar glitches, or potentially observable phenomena, such as gravita- 
tional waves from stellar oscillations or deformations of the star, involve incorporating global 
properties of the star, properties of the elastic crust, and the coupling between the two [lH6]. 

One goal of neutron star physics is to make use of the available observations to pro- 
vide data complimentary to those provided by terrestrial nuclear experiments in order to 
better understand the microphysics of dense matter |7j. In order to extract microphysical 
constraints from observations, the theoretical description of the observations should ideally 
use consistent models of the nuclear physics throughout the star in describing the prop- 
erties of the matter. Simple global properties of a neutron star such as its mass, radius 
and total moment of inertia are relatively insensitive to the properties of the crust itself, 
which contributes a small fraction of the total to each; in such cases the dominant theoret- 
ical ingredient is the nuclear equation of state (EoS) for the core [8|. However, in models 
where the crust-core interplay is important there will be a sensitivity to both global and 
microscopic core and crust properties. That the dominant uncertainties in the nuclear EoS 
can be usefully encoded in the nuclear symmetry energy is well known. The dependence of 
global stellar properties such as the radius, and crustal properties such as its thickness on 
the symmetry energy has also been well studied. However, little attention has been paid to 
the relationship between the symmetry energy, crust properties, and global properties in a 
consistent way. This paper aims to provide one step towards such a consistent approach. 

The outer and inner layers of the neutron star crust are defined by the absence and 
presence of a neutron fluid external to the nuclear clusters respectively. Additionally, a 
third layer is often predicted to be present at the bottom of the inner crust, consisting 
of exotic nuclear geometries collectively termed nuclear 'pasta'; this layer may be called 
the mantle [9]. The mantle can contain two distinct sub-layers: (1) the layer in which 
the pasta structures contain nuclear matter of non-zero proton fraction, surrounded by a 
fluid of pure neutron matter ('normal pasta'); (2) the inverse situation in which the pasta 
structures are bubbles containing pure neutron matter, surrounded by a higher density fluid 
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of nuclear matter with a non-zero proton fraction ('bubble pasta'). In the latter layer, 
protons themselves are macroscopically free. Whether these layers are to be included when 
modeling the bottom of the inner crust and the transition to the neutron star core will have 
important consequences for the physics of the crust-core interface. 

Global properties of neutron stars such as mass and radius are relatively insensitive to 
the crust EoS -P(p). However, there are a number of important crustal quantities that are 
sensitive to the crustal composition and thickness, both of which depend on the uniform 
nuclear matter EoS, the energy density of a homogeneous fluid of neutrons and protons as a 
function of density and proton fraction. In particular, as crustal matter comprises relatively 
low- or zero -proton fraction nuclear matter, crustal properties are sensitive to the symmetry 
energy of the nuclear EoS, which describes how the energy density changes as proton fraction 
varies [101 El- Experimental effort has focussed on constraining the magnitude of the 
symmetry energy at saturation density J, and the derivative of the symmetry energy with 
respect to density L. 

The density at which the transition from the crust to the core occurs has been shown to 
be sensitive to L [101 fT2] ? ], and hence the volume fraction occupied by the crust. The mass 
and moment of inertia fraction of the crust, important in the modeling of pulsar glitches [1], 
is sensitive to the transition pressure, which has been shown to correlate with the slope and 
curvature of the symmetry energy around saturation density [H]. The maximum and typical 
size of mountains on NS crusts [TSl US] and frequencies of torsional crustal oscillations [TT]- 
[20] depends on the thickness of the elastically rigid part of the crust and its shear modulus, 
which depends on the characteristic quantities which define the lattice: the inter-ion spacing 
a and nuclear charge Z [2j. The shear modulus throughout the inner crust [2TH23] is an 
important quantity in determining the frequencies at which the crust might undergo global 
oscillations. Various transport properties within the crust will also depend somewhat on the 
crustal composition through the densities of the various components, the electron fraction 
and the lattice spacing. The region of densities at which the neutrons are predicted to 
form a superfluid via the ^Sq channel has a model dependence; the possibility of layers of 
the crust in which the neutron fluid is normal exists, and will depend on the crust-core 
transition density, and the existence of a neutron superfluid in the pasta layers will also 
depend on their extent. Also, the melting temperature of the crustal lattice depends on the 
compositional parameters, and its magnitude relative to the superfluid transition density 
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will impact the dynamics of the crust. 

The 'pasta' phases [2H might have very different behaviors to the purely crystalline 
phase. The shear modulus of the pasta phases is still very uncertain, and will have a large 
effect on the crust-core boundary layer, of particular importance to the damping of r-mode 
oscillations [3] and the spectrum of torsional oscillations The phases of bubble nuclei 
may allow the direct URCA neutrino cooling process [9]. Thermodynamic and neutrino 
scattering properties of the pasta phases could be rather different |26H28] . The extent of the 
pasta phases in the crust will determine the importance of those phases to crustal processes 
and the dynamics of the crust-core coupling. The region of the inner crust occupied by the 
'pasta' phases is also strongly dependent on L, since, while the density of transition to pasta 
remains roughly constant in a given model pjOj, the crust-core transition density decreases 
with increasing L, reducing the density range of the pasta phases at higher L. 

The compressible liquid drop model (CLDM) |29] has proven to be a useful tool for the 
study of neutron star crust [HI |211 12SH2H1 EDHSH] • Like the liquid drop and droplet models 
of terrestrial nuclei, the CLDM focusses on the average, macroscopic properties of nuclear 
clusters in the crust such as mass and radius, whilst neglecting the quantum shell effects. 
The interior nuclear density is a free variable (hence compressible) rather than specified 
a priori as the saturation density of symmetric nuclear matter (SNM). An external gas 
of dripped neutrons, expected to coexist with the nuclear clusters in the inner crust, is 
calculated using the same bulk Hamiltonian as the matter inside the clusters. The neutron 
star inner crust is assumed to have the structure of a regular lattice, and the Wigner-Seitz 
(WS) approximation is employed in which the unit cell, which generally will have a cubic 
or more complicated structure, can be replaced by a spherical cell of the same volume. The 
W-S approximation is good so long as the nuclear size is much smaller than the cell size, 
but breaks down when the two sizes become comparable. Despite the two approximations, 
the CLDM provides a good description of the average microscopic properties of the neutron 
star crust, including the energetically preferred nuclear mass, size, WS cell size and density 
of dripped neutrons. Its amenability to quick calculation, and physical transparency, have 
made it the most widely used model for the computation of the EoS and composition of the 
inner neutron star crust and of inhomogeneous nuclear phases in core collapse supernovae 
HOI E]; it is therefore important to understand its limitations and accuracy given the 
experimental and theoretical uncertainty in its model parameters. 
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The purpose of this paper is to (i) provide a 'basehne' set of inner crust EoSs generated by 
varying the slope of the symmetry energy within hmits imposed by terrestrial experiments 
and theoretical pure neutron matter (PNM) calculations, together with their associated crust 
compositions, such that a suitable crust EoS can be found to match a core EoS with the same 
symmetry energy behavior, and (ii) systematically explore the remaining model dependence 
of the crustal composition, focussing on the experimental and theoretical uncertainties of 
the SNM EoS and the surface energy. In section II we outline the CLDM. In section III we 
describe our model for the bulk nuclear matter EoS, and the experimental and theoretical 
PNM constraints we will use. In section IV we discuss the range of parameters used in the 
proscription for the surface energy. In section V we present the set of baseline results, before 
probing the surface energy and SNM EoS dependence of the crust composition in sections 
VI and VII, and the dependence of the PNM EoS at intermediate densities in section VIII. 
In section IX we discuss how the EoS and surface uncertainties can impact important crustal 
properties such as its shear modulus and melting temperature, and we conclude in section 
X. 

II. THE COMPRESSIBLE LIQUID DROP MODEL 

We use the CLDM originally formulated by BBP [29j and updated to incorporate non- 
spherical nuclear shapes by Watanabe, lida and Sato [261 EH [39]. For completeness, we 
recall the main ingredients of the CLDM. Assuming a regular crystal structure for crustal 
matter, a repeating unit cell can be identified in which a single nucleus (or nuclear cluster) 
containing all the protons of the cell resides immersed in an external, uniform neutron fluid. 
A homogeneous gas of electrons is present to neutralize the positive protonic charge. We 
consider the five canonical nuclear geometries (the 'pasta' shapes): spherical, cylindrical, 
planar, cylindrical bubble and spherical bubble. The dimensionality of the shapes are spec- 
ified by a parameter d = 3,2,1 for spherical, cylindrical and planar geometries respectively. 
The Wigner-Seitz (WS) approximation is employed, in which the actual unit cell is replaced 
by an equally volumed unit cell of spherical, cylindrical or planar geometry depending on 
the pasta shape chosen. 

ece\lirc,X,n,nn) = w[nE(n,a;)+£exch + £^thick] +U5surf + ^i£Coul + (l-^^)'^n^(nn, 0)+5e(ne) (1) 
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E{n, x) is the energy per particle of uniform nuclear matter, ecoui the Coulomb energy 
including the proton-proton and proton-electron interactions, £surf is the surface energy and 
ee{ne) is the electron energy, u = {rf^/rcY is the volume fraction occupied by the nuclei or 
the bubble of radius tn, dimension d, and 

u nuclei , , 

(2) 

1 — u bubble 

is the volume fraction of matter of non-zero proton fraction, and distinguishes nuclei from 
bubble phases in equation ([T]). The independent variables are taken to be the WS cell 
radius r^, the density and proton fraction inside the nuclear cluster n,x and the neutron 
fluid density n^. The electron number density ng is related to these variables via the charge 
neutrality condition = unx. 

The global baryon number density is related to these variables through the relation 

= vn + {1 — v)n^ (3) 

from which we can define the fraction of nucleons in the neutron fluid = (1 — v)n^/ni, 
and in the nuclear cluster Xnuc = vn/ri]^ such that + Xnuc = 1- 

The Coulomb energy, which includes the lattice energy (the proton-electron contribution) 
is given in the WS approximation by 

e^c+L) = 2Tx{exnrNf fd{u) (4) 
1 



Uu) 



d + 2 



(5) 



_d-2 

The surface thickness correction from BBP [29] and the proton Coulomb exchange energy 
are given by 

4 

Sthickik.x) = --'Ke^w^x'^k^n (6) 

£exc/.(A:,a;) = -|^2i/3eV/3fcn (7) 

where k = (l.Svr^n)^/^ and w is a, distance representing the surface thickness, taken to be 
w ~ 0.75 fm. The energy density of the free electron gas is taken to be simply the free Fermi 
gas expression 



-hcn^/^xnu)^/^ (8) 



neglecting electron screening. 



The composition of the cell is obtained by minimizing the energy of the unit cell with 
respect to the free parameters. This produce four relations which correspond physically to 
mechanical, chemical and beta equilibrium of the cell plus the nuclear virial relation 

Ccoul 26surf ~l~ ^curv (9) 

which expresses equilibrium with respect to variation of relative volumes of the two phases 
in the unit cell. 

To simplify our model, we neglect electron screening and the effect of a neutron skin 
to the nuclear clusters. Shell effects both in the nuclei and the dripped neutron fluid are 
not taken into account in this classical model, and the WS approximation is expected to 
break down at the highest densities in the crust [:42|. In a later work, using a microscopic 
calculation, we will assess the effects of these approximations and missing physics using the 
results of this paper as a comparison. 



III. THE EQUATION OF STATE OF UNIFORM NUCLEAR MATTER 

Useful parameters characterizing the EoS of isospin asymmetric nuclear matter around 
SNM (proton fraction x = 0.5; 6 = 0) and the saturation density of SNM rig can be derived 
by expanding E{n,x) in a power series in the isospin asymmetry 6 = l — 2x and the density 
parameter x = ^^^^ 

E{n,x) = Eo{n) + S{n)6^ + ... (10) 
Eo{n) = Eo+'^Kox' + ... (11) 

S{n) = J + Lx+lK,yrr.X^ + ... (12) 

EpNM(n) ^ Eo{n) + Sin) (13) 

Eoin) = E{n, 0.5) is the binding energy per nucleon of SNM and S{n) = ^d'^E{n, x) / <95^^o 5 
is the nuclear symmetry energy. Kq is the incompressibility of SNM at saturation density. 
J = S^uq), L = dS{n)/dx\n=no and -ft'sym are the value of the symmetry energy, its slope 
and its curvature at saturation density. In particular, the pressure of pure neutron matter at 
sub-saturation densities, which plays a large role in determining the equilibrium composition 
of the crust, can be expressed as 

P=—[L + {Ko + K,,^)x + •••]• (14) 
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to the leading two orders. 

Throughout most of this paper we will mainly use the modified Skyrme-like (MSL) pa- 
rameterization of the nuclear matter EoS E{n, 6) (see appendix A) as our description 
of uniform nuclear matter as a function of density and isospin asymmetry. Like the Skyrme 
parameterization, the MSL model has a number of free parameters; these parameters can 
be analytically related to the properties of uniform nuclear matter at saturation density, 
allowing a smooth variation of, e.g., the symmetry energy at saturation J and its slope L. 
For comparison, we will also use a similar phenomenological EoS whose form was originally 
written down by Bludman and Dover [44] (which we will refer to as BD, see also appendix 
A), which was later modified and used to study finite nuclei and inner crust composition by 
Oyamatsu and lida [101 HH], and a selection of Skyrme EoSs [46j whose basic properties are 
given in Table 1. 

We note that for all of these EoSs, the symmetry energy is calculated in full without 
the use of the parabolic approximation (PA) in which one truncates the expansion in equa- 



tion (10) at second order, thus obtaining the symmetry energy as S{n) = Epj^Mi^) — Eq^u). 
Although, in many circumstances, the parabolic approximation is good enough, since the 
coefficients of fourth order and higher are small, it has been shown ^I3l [H] that for se- 
lected crustal properties, most notably the crust-core transition density, the PA can lead to 
divergent predictions compared to the full EoS. 



TABLE I: Skyrme parameters. J, L, and Kq are in MeV, uq in fm 





J 


L 


Ko 


no 


SLy4 


32.04 


46.03 


230.31 


0.16 


SkM* 


30.06 


45.79 


216.97 


0.161 


BSkl4 


30.03 


43.94 


239.74 


0.159 


BSkl6 


30.02 


34.89 


242.08 


0.159 


NRAPR 


32.82 


59.69 


226.06 


0.161 


SII 


34.15 


50.07 


341.19 


0.148 


SkM 


29.54 


60.51 


248.41 


0.16 


SkX 


31.11 


33.14 


271.48 


0.156 
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A. Experimental constraints on the symmetry energy 

The magnitude of the symmetry energy at saturation density J is constrained mainly from 
nuclear mass model fits to experimental data . Care must be taken when interpreting 

predictions for J, since often different definitions for the symmetry energy at saturation are 
used, for example differing in the use of the parabolic approximation. Moreover, in mass 
model fits, J depends significantly on what surface symmetry energy is used, e.g. [51]. We 
take conservative ranee 25 < J < 35 MeV. 

Constraints on the density dependence of the nuclear symmetry energy have been ob- 
tained in analyses of a variety of nuclear experimental data E3]- Recent modeling of 
isospin diffusion in heavy ion collisions involving ^^^Sn and ^^^Sn extracted constraints of 
62< L <107 MeV using the IBUU04 transport model [Ml |55] and 45< L <103 MeV using 
the ImQMD molecular dynamics model [53] respectively. The extraction of the range of 
L is model dependent; the two ranges quoted here, whereas they overlap significantly, are 
from the two different transport model analyses. A study of isoscaling data from multi- 
fragmentation reactions yields L ~ 66 MeV [SB]- A constraint from the analysis of pygmy 
dipole resonance (PDR) data gives 27< L <60 MeV [57j and from the analysis of the sur- 
face symmetry energies of nuclei over a wide range of masses gives 75 < L <115 MeV [58] . 
A measurement of neutron skins of a wide mass range of nuclei has led to an estimate of 
25< L <100 MeV [32] • Combining data on the neutron skin thickness of Sn isotopes and 
isospsin diffusion and double n/p ratios in heavy ion collisions, a range of 40 < L < 70 MeV 
is obtained [60], while constraints on the global nucleon optical potential from nucleon- 
nucleus reactions and single particle energy levels of bound nuclear states leads to a range 
30.2 < L < 75.2 MeV [61J. For the latest and more extended discusssions on the experi- 
mental constraints on J and L, we refer the reader to the NuSYMll (Nuclear Symmetry 
Energy Workshop 2011) website |62]. We take as a conservative range 25 < L < 115 MeV, 
noting however that more recent studies favor the lower half of this range [SHI E3| • 

As a representation of the above uncertainties, we focus on the PNM EoS as this enables 
us to make contact with microscopic calculations the predictions of which, at low densities, 
are now quite robust. For a given model, the SNM EoS {5 = 0) and the symmetry energy 
fixes the PNM EoS. The PNM EoS for the MSL and BD models at a constant value of 
J = 35 MeV with L = 25, 70, 115 MeV are plotted in the left panels of Fig. 1. For constant 
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FIG. 1: (Color online) PNM energy per particle versus baryon number density for the MSL and BD 
EoSs with J=35 MeV (left), constraining the MSL and BD EoSs to fit the low-density end of the 
PNM EoS from microscopic calculations (middle) and for a selection of Skyrme EoSs (right). The 
gray and dotted bands correspond to the microscopic calculations of Gandolfi et al [69] (GCR) and 
Hebeler and Schwenk |68j (HS) taking into account uncertainties in the three-nucleon interaction, 
while the boxed region at low density comes from the calculations of Schwenk and Pethick |66| 
(SP). 



J, the MSL and BD EoSs match closely for low densities and high values of L. 



B. Theoretical constraints on the pure neutron matter EoS 

Theoretically, the PNM EoS has recently been well constrained at low densities through 
quantum Monte-Carlo, Green's function Monte-Carlo, chiral effective field theory and vari- 
ational chain summation techniques [631 - I68] . These constraints, which are expected to be 
quite robust, may be taken as a constraint for our uniform nuclear matter EoS. We show the 
PNM constraint from Schwenk and Pethick (SP) [66] as the red box at low densities; results 
since have converged on the lower bound of that box. Calculations at higher sub-saturation 
densities introduce uncertainties which are dominated by the theoretical uncertainties in 
the three-nucleon force; two recent calculations with estimated theoretical error bars have 
been performed [Ml [69] (labelled HS, GCR hereafter, respectively), and are indicated by the 
shaded bands in Fig. 1. 

In the left panel of Fig. 1, one can see that our phenomenological EoSs at a constant 
value of J are consistent only for a very limited range of L (e.g. L ^ 70 MeV for J = 35 
MeV for the MSL EoS). We thus explore the effect of constraining our phenomenological 
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FIG. 2: (Color online) Correlation between J and L obtained by requiring the MSL EoS pass 
through the upper bound from low density PNM calculations (upwards pointing triangles) and the 
lower bound (downwards pointing triangles); the gray shaded region indicates the region consistent 
with the PNM constraints. The yellow band is obtained by applying the Hugneholtz-van-Hove 
theorem to optical potential data [70], the black lines from PNM calculations at higher densities 
(red plusses - chiral perturbation theory [6&]; red crosses - QMC methods [6^) and the blue lines 
from constraints on the saturation properties of asymmetric matter from mass fits [lOj . 

EoSs using PNM calculations as follows. 

We constrain our MSL and BD EoSs given a reference low density PNM EoS for a given 
value of L by varying the magnitude of the symmetry energy at saturation J for that L until 
our phenomenological PNM EoS passes through the reference EoS. For another value of L, J 
must be readjusted so that we match the reference EoS again. In the middle panel of Fig. 1, 
we plot the MSL and BD PNM EoS following such procedure, with the reference EoS at low 
densities taken to be the results of HS and GCR which closely coincide with the lower bound 
of the SP box. We show results for L = 25, 70 and 115 MeV. For a given value of L, and 
applying the low density PNM constraint, the MSL and BD EoSs differ significantly, with 
the BD PNM energy being appreciably lower for most of the density range up to saturation 
(corresponding to a lower symmetry energy J); thus the BD EoS will give a lower pressure 
of PNM. Finally, we display a selection of Skyrme PNM EoSs in the right panel of Fig. 1; 
the selection was made based on the criterion that they pass through the SP constraint. 

One trivial consequence of applying this constraint is that we are no longer free to vary J 
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independently of L; we select only those EoSs in J-L space which conform to the low density 
PNM constraints, thus imposing a correlation between the symmetry energy and its slope 
at saturation density. This correlation is shown in Fig. 2 for the MSL EoS by the upright 
black triangles (for the upper bound of the SP constraint) and the lower black triangles (for 
the lower bound), and for the BD EoS by the squares and diamonds respectivley. Taking 
the upper or lower bounds of the SP constraint produces a different offset to the correlation 
without changing its slope. The existence of such a correlation is in part due to the variation 
of the PNM EoS with density; all theoretically derived energy- density functions are simple 
enough to admit such a correlation. For example, a correlation derived from the low density 
PNM constraint might not exist if the energy of PNM is not a monotonic function of density 
in the region up to saturation density. This, however, would imply a region of density in 
which PNM became stable, a phenomenon that we would expect to have seen indications 
of in nuclear experiments. Given a smooth, monotonically increasing function of density, 
a constraint at the low density end will give a correlation between its value at the higher 
density region, and its slope there. 

Although our correlation is simply a result of applying the PNM constraint, such a corre- 
lation has been noted in several other works, also plotted in Fig. 2. The PNM calculations 
can be extended towards saturation density and beyond with the inclusion of uncertain 
three-body forces; given an empirical value for the symmetry energy at saturation, one then 
obtains a value for L. Correlations in the J — L plane are shown for two such recent studies 
EH]- The Hugenholtz- Van-Hove theorem predicts a relation between J and L whose 
uncertainty, in the symmetry potential, can be related to nucleon optical potentials; in the 
analysis of Xu et al [70j , the yellow shaded region in Fig. 2 is predicted. Oyamatsu and lida 
[To] , hereafter 01, obtained a constraint from fits to nuclear masses and charge radii using 
the BD EoS of 5 = 28 + 0.075 L. Finally, our selection of Skyrmes are also plotted. As can 
be seen, there is some variation in the slope of the correlation and its offset from the L axis, 
but the studies so far predict correlations that favor the low-L, high J half of the plane, in 
which our correlation also falls. The 01 results show that two different correlations can be 
obtained for the same EoS depending on the data used to obtain it. 
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IV. THE SURFACE ENERGY 



The surface energy density can be written in terms of the surface tension cJs as 

£surf — ) l-'-'-'J 

rN 

where tn is the radius of the nuclear cluster (half-width in the case of planar clusters) and 
d the dimensionality of the cluster. In order to take into account the effect of the curvature 
of the nuclear surface, a curvature correction is added: 

(Ts CTs + (c? - l)(Tc/rN, (16) 

where cXc is the curvature tension. 

Analogously to nuclear matter, the surface and curvature tension may also be expanded 
about its isospin symmetric value: 

cTsix) = ao + as5^, (17) 

cTc(x) = cTco + crc,sS^, (18) 

where as = ^d'^(Ts/d6^\x=o.5 is the symmetry surface tension; the symmetry curvature tension 
is defined similarly. 

The first formulations of the CLDM parameterized the surface energy in terms of the 
proton fraction x, the bulk nuclear and neutron matter densities n, n^, and the difference in 
energy per particle between the two bulk phases P^HSTl inspired by one- and two-fluid 
Thomas- Fermi theory. It was noted that, in fact, the surface energy is to be evaluated along 
the curve of coexistence between the bulk nuclear and pure neutron matter [7H [72] and 
therefore the quantities upon which it depends not all independent variables. 

At zero temperature, the surface tension may be parameterized in terms of one variable, 
which is most conveniently chosen to be the proton fraction x. The inclusion of a neutron 
skin accounts for the dependence of surface properties on the bulk nuclear densities. In this 
work we omit the neutron skin; the effects of the skin will be examined in a subsequent 
study. Here we use the parameterization first used by Lattimer, Pethick, Ravenhall and 
Lamb (LPRL) |7T1 172], which was subsequently generalized by Lorenz and Pethick |73l [7i] : 

2^+1 + b , , 

" {l-x)P 
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FIG. 3: Surface symmetry parameter cJsym versus symmetry energy J for four values of the surface 
symmetry - volume symmetry correlation. 



from which we can derive 



The curvature tension takes the form [731 El] 

fTc = fT,^a(/3-x)+7a;^ (21) 

The parameters ao, o-q^c, a, /3, 7 are adjusted to fit calculations of the interface between 
two different phases of semi-infinite nuclear matter using the Thomas-Fermi and Hartree- 
Fock methods 



A. Constraints on the surface energy 

Experimentally, o"o, o"c and are determined largely by fitting liquid drop models to 
nuclear mass data [HI HSJ ED I75H75] . The values obtained are rather dependent on the 
formulation of the liquid drop model. The strength of the surface and curvature tensions 
for SNM are relatively well determined to be in the range of 1.0 < ctq ^ 1-1 MeV fm~^, 
0-5 ^ o"c,o ^ 0.6 MeV fm~^; in our model we fix (Tq = 1-1 MeV fm~^ and (7c,o = 0.6 MeV 
fm~^. The surface symmetry energy is rather more model dependent; that uncertainty can 
be systematically encoded in the correlation which exists between the surface symmetry 
energy and the bulk symmetry energy J. The correlation emerges both from fits to nuclear 
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FIG. 4: (Color online) Form of the surface and curvature potentials used in this work for 3 different 
values of the surface symmetry energy and 2 different parameterizations of the low proton fraction 
behavior. 



masses and from calculations of the properties of plane interfaces; the slope of the correlation 
depends on the model used [771 EH] • We incorporate this correlation as follows: 

E n'^ 

08 = (^^^yi^[(0-046c + 0.01)Esym - c] (22) 

where c gives the slope of the correlation. In Steiner et al [77] a comprehensive study was 
conducted of the extraction of the surface symmetry energy from mass model fits and micro- 
scopic calculations of plane interfaces between two different phases of nuclear matter. The 
surface symmetry energy was extracted using two different, thermodynamically consistent 
formulations of the surface thermodynamic potential in the mass model; the so-called '/in' 
and '/Iq,' approaches which differ in the way the neutron skin is taken into account. The best 
fits obtained for the '/in' and '/i^' approaches are given by c = 3.14 and c = 4.8 respectively 



in equation (22); the full region of the surface symmetry - bulk symmetry energy plane can 



be spanned by c = 2.0 to c = 7.0. In our model, we impose the correlation in equation (22) 
and take for our baseline results c = 4.8. We can then explore the effect of the correlation 
systematically by varying c over the full range from c = 2.0 to c = 7.0. In Fig. 3 this range is 
depicted. In Fig. 4 we show the form for the surface and curvature potentials for 3 different 
values of the surface symmetry energy as corresponding to J = 30 MeV, c = 2, 4.8 and 7 
by the black solid, long-dashed and short dashed lines respectively. 

Of particular importance is the behavior of the surface energy at low proton fractions. 



determined in equation (19) by the parameter p. Although the fits to plane interface calcu- 
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lations suggest a value of p = 3, these were performed using Skyrme interactions. Recent 
work on neutron drops suggests that the gradient terms in the Skyrme interactions underes- 
timate the surface energy for pure neutron matter [TH] . We will therefore also vary p in the 
range 2 - 4 to determine the relative importance of the low proton fraction surface energy. 
In Fig. 4 we also show the form of the surface and symmetry energy for p=2 and 4 for as 
= 3.6 corresponding to c = 4.8, J = 30 MeV. The middle panel shows a close-up of the 
low proton fraction region of the surface potential, emphasizing the fact that p = 4 gives a 
higher surface tension in this region, p = 2 a lower surface potential. 

V. BASELINE RESULTS 

Our baseline SNM EoS and surface parameters are summarized in Table 2. We then vary 
the slope of the symmetry energy L from 25-115 MeV under two different constraints: 

• Constant symmetry energy at saturation density J. Varying J over the range 25 - 35 
MeV generates a new sequence of EoSs over the range of L; we will use as a shorthand 
label for such a sequence at, say J = 35 MeV, 'J35' 

• Demanding a specific PNM EoS at low densities; then J will vary with L according 
to the correlations given in Fig. 2. We will concentrate on the lower bound of the SP 
constraint outlined in Figs. 1 and 2; we will label the sequence of crust EoSs generated 
by such a sequence as 'PNM'. 

In our model, the region of the J, L parameter space consistent with both the conservative 
experimental range of J = 25 — 35 MeV and the theoretical PNM constraint corresponds to 
the shaded gray region in Fig. 2; the corresponding region for the various crustal properties 
will also be shown as a shaded region in subsequent plots. We also note that, since the recent 
results from PNM calculations are expected to be quite robust, and give a low density PNM 

TABLE II: Baseline parameters 
Kq (MeV) no(fm-3) ao (MeV fm-^) crc,o(MeV fm"!) c p 
240 0.16 1.1 0.6 4.8 3 
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FIG. 5: Crust-core and spherical nuclei-pasta transition densities (left) and pressures (right) versus 
L. Results are displayed for the 3 constant-J sequences 'J35', 'J30' and J25' and the two sequences 
constrained by the low density PNM EoS 'PNM'. The shaded region correspond to the region 
consistent with the PNM constraint and 25 < J < 35 MeV. 



EoS at the lower bound of the SP constraint, the 'PNM' sequence gives the most stringent 
constraint on crust properties within our model. 

We will concentrate on the crust-core, and 'normal' (spherical) nuclei - pasta (cylindrical 
phase) transition densities and pressures (ncc, ^p. Pec, -Pp respectively), and on the variation 
throughout the inner crust of the pressure P and the following four compositional parameter: 
with nn the local free neutron density and the overall baryon density of matter, 

• V, the volume fraction of clustered matter (the component which contains the protons), 

• Xn, the density fraction of free neutrons, = (1 — v)n.a_/n\,, 

• rc, the radius of the WS cell, and 

• x, the local proton fraction of clustered matter. 

Any other compositional quantities can be expressed in terms of these four; for example, 
the mass and charge number of nuclei A = 4'7rr^(l — X^)/3, Z = xA over the range where 
spherical nuclei exist (it is not defined for the pasta phases as the nuclear clusters are not 
localized to one WS cell). 

In Fig. 5 we show the variation of the crust-core and spherical nuclei-pasta transition 
densities and pressures over the range of L. Results are plotted for EoSs using the PNM 
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low constraint and also constant J = 25, 30, 35 MeV. For the crust-core transition we also 
show results for the PNM high constraint for reference. 

Whether we choose the PNM constraint or constant J constraint, the crust-core transition 
density shows the expected negative correlation with L, although for the PNM constraint 
the correlation is much less pronounced; for PNM low, rice ~ 0.07 — 0.095 fm~^, a variation 
of 0.25 fm-2; for constant J = 35 MeV rice ~ 0.05 - 0.012 fni'^, a variation of 0.07 fni'^. 
The crust-core transition density also clearly depends on J with a positive correlation, and a 
variation of ~ 0.3 fm~^ from J = 25 — 35 MeV independent of L; a higher symmetry energy 
favors a larger proton fraction in the nuclear clusters, lowering their bulk binding energy to 
a greater extent than the increase in symmetry energy, thereby making the clustered matter 
energetically favorable to higher densities. For the 'PNM' sequence, J increases with L, 
and so the decrease of the transition density with L is compensated slightly by the increase 
of the transition density with J, making the rice — L correlation shallower. Indeed, with a 
sufficiently steep correlation between J and L, the correlation can vanish or even be reversed; 
if we were to take the HVH correlation of Fig. 2, we would see no correlation between ^cc, 
and L. The shaded region, consistent with experiment and the SP constraint, encompasses 
crust-core transition densities of 0.08 — 0.12 fm~^. 

Relative to the crust-core transition density, the spherical nuclei-pasta transition density 
shows very little variation with L and J, tracing out a thin band rip ~ 0.05 — 0.06 fm~^. 

The transition pressures show very different correlations with L depending on whether 
we choose the PNM constraint (positive correlation) or the constant J constraint (negative 
correlation). However, it can be seen that transition pressure has a positive correlation 
with J; for the PNM constraint J increases with L, and so overall the correlation is deter- 
mined by the competition between increasing Pec with J and decreasing Pec with L. The 
spherical-pasta transition pressure Pp also increases with L for the PNM constraint, whereas 
the variation with L at constant J varies more weakly with L with the same dependence, 
phenomenologically, as the transition density. 

The crust-core and spherical-pasta transition pressures consistent with PNM low and 
25 < J < 35 MeV both vary by more than a factor of 2, Pec ~ 0.35 - 0.85 MeV fm~^, 
Pp fti 0.15 - 0.3 MeV fm'^. 

Note that the shaded region is bounded by the 'PNM' and 'J35' sequences for L < 70 
MeV; these two sequences will be taken as the bounding regions in what follows. 
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FIG. 6: (Color online) Pressure of matter scaled by n^^^ versus baryon number density in the inner 
crust. Results are displayed for the L = 25,70 and 115 MeV members of the 'J35' and 'PNM' 
sequences. The left panel displays the results over the whole density range of the inner crust, 
while the right panel shows the lower density region up to the transition to pasta. The loci of the 
transition densities over the range L = 25 — 115 MeV are displayed for the 'PNMlow' sequence 
(larger black circles) and the 'J35' sequence (smaller blue circles). The shaded region has the same 
meaning as in Fig. 5. 



We now look at the ranges of compositional parameters predicted in the basehne model. 
Results are shown for the 'PNM' and 'J35' sequences for L = 25,70 and 115 MeV. Note 
that the two sequences coincide at L = 70 MeV. The left panels show the full density range 
on a linear scale, which emphasizes the values within the range of the pasta phases. The 
loci of the crust-core transition densities are shown for both PNM low (larger black circles) 
and constant J (smaller blue circles). The right panels show the range of densities up to the 
transition to pasta; when the predictions using the different L and J values differ significantly 
at the lower densities, a logarithmic scale is used to clearly display those differences. The 
loci of the spherical-pasta transition densities are shown for both PNM low (larger black 
circles) and constant J (smaller blue circles). 

Fig. 6 shows the total pressure of matter throughout the crust, scaled to n^^^; a free, 
relativistic neutron gas would result in a constant. Most models, throughout the middle 
regions of the inner crust, predict a roughly constant trend too. The main deviation comes 
from J = 35 MeV, L = 115 MeV, values which predict the slope for PNM energy per 
particle versus density, and hence the pressure of PNM, to approach zero at low densities. 
As expected, since the pressure comes to be dominated by the free neutrons in the crust, a 
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FIG. 7: (Color online) The volume fraction v occupied by nuclei/nuclear clusters (top), and density 
fraction occupied by free neutrons (bottom), versus baryon number density in the crust. The 
same selection of EoSs from Fig. 6 is used, and the transition loci are displayed in the same way. 
The shaded region has the same meaning as in Fig. 5. The left panels display the results over the 
whole density range of the inner crust, while the right panels show the lower density region up to 
the transition to pasta. 



higher value of L leads to a higher pressure throughout most of the crust. 

Fig. 7 shows the volume fraction of clustered matter and the density fraction of free 
neutrons. The first thing to note is that the volume fraction reaches > 0.8 at the crust-core 
transition for the shaded region; this emphasizes the fact that the WS approximation is 
certainly not valid at the highest crustal densities. For the PNM constraint, the crust-core 
transition fraction weakly depends on L, ranging from 0.8 (L=115 MeV) to close to 1.0 
(L=25 MeV). For J = 35 MeV, the crust-core transition fraction varies much more widely, 
from close to 1 (L=25 MeV) down to 0.1 (L=115 MeV). The spherical-pasta transition 
fractions vary over a much smaller range from ^ 0.11 - 0.14 for PNM low and ~ 0.08 - 0.12 for 
J = 35 MeV. A useful comparison is with the simple estimate from the Bohr- Wheeler fission 
instability criterion (see, e.g. [BO]) of 0.125, a value which closely matches the range predicted 
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FIG. 8: (Color online) Local proton fraction of clusters x (top) and WS cell size rc (bottom) 
versus baryon number density in the crust. The same selection of EoSs from Fig. 6 is used, and 
the transition loci are displayed in the same way. The shaded region has the same meaning as in 
Fig. 5. The left panels display the results over the whole density range of the inner crust, while 
the right panels show the lower density region up to the transition to pasta. 



by our baseline model. This value is indicated on Fig. 7 by the horizontal line. Below that 
density, the fraction varies by up to 0.05 from L=115 MeV (lower values) to L = 25 MeV 
(higher fractions), an easily understood correlation: higher L implies higher neutron matter 
pressure at saturation density and below, resulting in the pressure equilibrium condition 
favoring smaller nuclei. 

Mirroring the volume fraction at the crust-core transition, varies from close to 0.0 up 
to 0.2. The fraction rapidly converges to zero at low densities as the density of free neutrons 
vanishes. Xn peaks around 0.01 - 0.03 fm~'^ with values between 0.8 - 9.6, with higher L 
being correlated higher Xn. The transition to pasta occurs close to the peak value of Xn for 
all L and J. 

Fig. 8 shows the cell size rc and the proton fraction of clustered matter x. rc varies 
smoothly up to the transition to pasta, and then proceeds through a series of discontinuous 
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jumps as matter transitions through the various nuclear shapes. In practice, these jumps may 
be smoothed out by the existence of intermediate shapes not considered in this work |HT] . 
The negative correlation with L is again explained by the pressure equilibrium condition. 
In the right panel, it can be seen that the variation of rc over the density range up to the 
transition to pasta consistent with our PNM and constant J constraints is constant at ~ 8 
fm. 

X generally decreases with increasing density over the whole density range. For J = 35 
MeV, the crust-core transition fraction varies by about 0.03; from L = 25 MeV, Xcc decreases 
from 0.06 down to 0.03 at L = 70 MeV, then increasing back up to 0.06 at L = 115 MeV; 
for PNM the variation is similar, starting at Xcc=0.025 for L = 25 MeV and increasing 
with L up to 0.06 at L = 115 MeV. In the lower density region the variation in x remains 
approximately constant at around 0.05. At a given density, higher J and higher L correlate 
with higher x. Higher J favors a higher proton fraction; higher L favors smaller, denser 
nuclei and nuclear clusters; as the symmetry energy increases with density, denser nuclei 
will also tend to favor higher x. 

VI. SURFACE ENERGY DEPENDENCE 

Having established how the crustal parameters behave in our baseline sequences, we begin 
to examine how those results depend on the various other model parameters, starting with 
the surface and curvature parameters. In the following figures, the baseline results will be 
represented as the shaded regions, bounded by the 'PNM' and 'J35' sequences and L < 70 
MeV. 

In this section we consider the same set of transition and compositional parameters as 
in the previous section. We examine changes in the behavior of the surface and curvature 
energies at proton fractions near 0.5 (as controlled by the surface symmetry parameter as 
through the correlation with the bulk symmetry energy J, controlled by the parameter c) 
and low proton fraction, as controlled by the parameter p. We vary c from 2.0 to 7.0 (which, 
for any given value of J corresponds to an increase in as by a factor of 3.5); note, a higher 
value of c, and therefore cr^, gives a lower surface energy at a given proton fraction X. 
We also examine the cases of p = 2 and p = 4 as mentioned in section IV. A. p = 2 gives 
a substantially higher surface energy at low proton fractions, while p = 4 is substantially 
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FIG. 9: Crust-core and spherical nuclei-pasta transition densities (left) and pressures (right) versus 
L for the 'PNM' (thick lines) and 'J35' (thin lines) sequences. The top two plots show the results 
with c = 2.0 (solid lines) and c = 7.0 (dashed lines); the bottom two plots show the results with 
p = 2 (solid lines) and p = 4 (dashed lines). The baseline results are shown for reference as the 
shaded regions. 



lower (see Fig. 4). 

Fig. 9 displays the transition densities and pressures. The top two figures show the 'PNM' 
and 'J35' sequences with c = 2.0 and c = 7.0; the lower two figures show the two sequences 
for p = 2.0 and p = 4.0. A higher surface and curvature energy will tend to favor a lower 
transition density, as the energy cost of having a surface is higher; for a given value of L 
on any of the sequences, we would thus also expect a higher transition pressure. This trend 
should be the opposite for the spherical nuclei-pasta transition; the pasta shapes generally 
have a larger surface, so a higher surface energy would favor the spherical nuclei to higher 
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FIG. 10: (Color online) Pressure of matter scaled by n^/^ versus baryon number density in the 
inner crust. Results are displayed for the L = 25 and 70 MeV members of the 'J35' and 'PNM' 
sequences. The left panel displays the results at c = 2 and c = 7; the right panel for p = 2 and 
p = 4. The baseline region is shown for reference. 



densities. Thus a higher surface energy decreases the density range over which pasta nuclei 
are energetically favorable. 

This is indeed what we see; the c = 2.0, p = 2 sequences have lower crust-core transition 
densities and higher spherical-pasta transition densities compared to the baseline (c = 4.8) 
results, while the c = 7.0, p = 4 sequences have elevated crust-core transition densities, and 
lower spherical-pasta transition densities. 

The effect on the transition densities is relatively small when c is varied, amounting to 
~ 0.005 fm~^ around the baseline results for the most part. The variation with p is more 
dramatic; for p = 4 the difference with the baseline results is still relatively small (< 0.005 
fm~^) but for p = 2 the crust-core transition density is significantly lower then the baseline 
results, falling by 0.02 fm~^ throughout the 'PNM' sequence, and by up to that amount 
for the 'J35' sequence at high L. The effect on the spherical nuclei-pasta transition is not 
as dramatic; the end result is that the region in which the pasta shapes are present is 
significantly reduced; the corresponding effect on macroscopic crust properties is discussed 
in section IX. 

For the compositional parameters throughout the inner crust, we display results for the 
L = 25 MeV and L = 70 MeV members of the 'J35' and 'PNM' sequences, which together 
bound the shaded regions for the baseline results. Fig. 10 shows the pressure scaled by n^/^. 
Fig. 11 shows the volume fraction v and density fraction occupied by free neutrons, and 
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FIG. 11: (Color online) The volume fraction v occupied by nuclei/nuclear clusters (top), and density 
fraction occupied by free neutrons (bottom), versus baryon number density in the crust. Results 
are displayed for the L = 25 and 70 MeV members of the 'J35' and 'PNM' sequences. The left 
panel displays the results at c = 2 and c = 7; the right panel for p = 2 and p = 4. The baseline 
region is shown for reference. 



Fig. 12 the local proton fraction of clusters x (top) and WS cell size vq] in each case the 
baseline results (shaded regions) are compared with the upper and lower limits of the surface - 
bulk symmetry correlation c = 2.0, 7.0 (left hand panels) and the low-proton fraction surface 
parameter p = 2,4. The pressure exhibits very little dependence on either of the surface 
parameters, dominated as it is by the pressure of the dripped neutrons. The volume fraction 
shows little dependence on c throughout most of the crust; at higher densities it reaches the 
same maximum values as the baseline at a slightly lower (higher) density for c=2.0 (7.0). 
Throughout most of the crust, p has little effect on the volume fraction; however, just as 
p=2 gives a large reduction in the transition density and pressure, the volume fraction is 
also reduced significantly, down to about 0.3, reflecting the fact that the increased surface 
energy becomes prohibitive when the clusters occupy a much smaller fraction of the space. 
The free neutron fraction exhibits the same dependence on c and p, with a significant effect 
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FIG. 12: (Color online) Local proton fraction of clusters x (top) and WS cell size rc (bottom) 
versus baryon number density in the crust. Results are displayed for the L = 25 and 70 MeV 
members of the 'J35' and 'PNM' sequences. The left panel displays the results at c = 2 and c = 7; 
the right panel for p = 2 and p = 4. The baseline region is shown for reference. 



only for p=2. 

Similarly, c and p have relatively little effect on the proton fraction throughout most of 
the crust; a lower value of c or p (higher surface energy at high and low proton fractions) is 
associated with a marginally higher proton fraction at lower densities, as an increased proton 
fraction lowers the energy of the bulk matter inside the nuclei (makes it more bound) to 
compensate for the increased surface energy. In contrast, the size of the WS cell does exhibit 
a significant dependence on the surface parameters, shifting the baseline sizes on the order 
of 5 fm higher (lower) for lower (higher) values of c and p throughout most of the inner 
crust. Since the volume fractions are relatively unchanged, this means that the nuclear sizes 
are also higher (lower). This can be understood from the competition between the surface 
and bulk energies: raising the surface energy (by lowering c or p) can be compensated by 
increasing the contribution to the total energy density of the bulk phases of matter relative 
to the surface contribution. 
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FIG. 13: Crust-core and spherical nuclei-pasta transition densities (left) and pressures (right) 
versus L for the 'PNM' (thick lines) and 'J35' (thin lines) sequences. The top plots show the 
results for the saturation densities no = 0.14 fm^^ (solid lines) and uq = 0.17 fm~^ (dashed lines), 
with the densities and pressures scaled to no; the bottom two plots show the results with Kq = 220 
MeV (solid lines) and Kq = 260 MeV (dashed lines). The baseline results are shown for reference 
as the shaded regions. 



VII. DEPENDENCE ON INCOMPRESSIBILITY AND SATURATION DENSITY 
OF NUCLEAR MATTER 

We now examine the impact on the transition properties and composition of varying 
the properties of SNM at saturation densities via the incompressibility Kq and saturation 
density no within experimentally constrained values, relative to our baseline results. The 
ranges of variation used are 220< Kq <260 MeV and 0.14< no <0.17 fm~^. 
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FIG. 14: Pressure of matter scaled by n^^^ versus baryon number density in the inner crust. Results 
are displayed for the L = 25 and 70 MeV members of the 'J35' and 'PNM' sequences. We show 
results for uq = 0.14(0.17) fm~^ (thick, (thin) lines) versus density scaled to uq in left panels and 
Kq = 220(260) MeV (thick (thin) lines) in the right panels. The baseline region is shown for 
reference. 




FIG. 15: The volume fraction v occupied by nuclei/nuclear clusters (top), and density fraction 
occupied by free neutrons (bottom), versus baryon number density in the crust, no(ii^o) a-^e in the 
left (right) panels in the same way as Fig. 14. 
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FIG. 16: Local proton fraction of clusters x (top) and WS cell size rc (bottom) versus baryon 
number density in the crust. The left and right panels show the results for the variation in uq and 
Kq in the same way as Fig. 14. 



Decreasing (increasing) Kq and no results in a decrease (increase) in the crust-core tran- 
sition density and pressure. In Fig. 13 we plot the transition densities (left) and pressures 
(right) versus L for % = 0.14,0.17 fm"^ (top two plots) and Kq = 220,260 MeV (bottom 
two plots), for the 'PNM' and 'J35' sequences. Results are compared with the baseline 
results (shaded area). 

Varying no simply rescales the EoS; we thus plot the transition densities and pressures 
scaled to the new saturation density in the top two plots. Plotted in this way, Ucc and 
Pec show little variation from the baseline results, indicating that, to good accuracy, the 
saturation density simply rescales them linearly. The lower (higher) Kq values give lower 
(higher) values for n^c and by up to ~ 0.005 fm~^; for the transition to pasta the effect 
is greatest at high L. Similarly, the pressures are lowered (raised) for lower (higher) Kq; the 
effect is greater at high L, up to ~ 0.005 fm~'^ for Pec, the effect on the pasta transition is 
small. 

In Figs. 14-16 we show the effect from the variation of Kq and no on the pressure and 
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composition throughout the crust. Again, when examining the effect of the saturation 
density, we scale the crust density by uq. P/n^^^ shows httle variation with no (once the 



density is rescaled) and Kq, although equation 13 shows a first-order dependence on Kq, 
its variation is much smaller than that of L. v and show a very small dependence for 
the same reason, as they are largely determined by the symmetry energy and pressure. 
Similarly x, determined in large part by the symmetry energy at a given density, shows 
negligible dependence. The WS cell size tq is the only parameter to show a significant 
dependence, increased by ~ 5 fm relative to the baseline results by lowering the saturation 
density. 



VIII. DEPENDENCE ON EOS 



So far we have examined the dependence of crustal composition on the properties of the 
EoS at saturation density L, Kq and uq while fixing the absolute value of the energy of PNM 
either at saturation density (by keeping J constant) or at very low densities (the PNM 
constraint). However, as Fig. 1 shows, this still allows some freedom of the PNM EoS at 
intermediate densities; for example, the MSL and BD EoSs, constrained by the same low 
density EoS, differ significantly at intermediate densities, the BD EoS giving a lower energy 
of PNM (and therefore pressure) for a given value of L, a fact that shows up in the difference 
in the correlations obtained in Fig. 2. We now examine the effect of these differences. 

In Fig. 17 the crust-core transition densities and pressures are shown for the 'PNM' and 
'J35' sequences of the BD EoS, as well as for the Skyrmes listed in Table 1, compared to 
our baseline MSL results. For the BD EoSs, the baseline quantities in Table 2 are used; the 
Skyrmes use the same surface parameters, and their saturation EoS properties are given in 
Table 1. There is little difference between the MSL and BD baseline results for the transition 
densities, and the Skyrmes fall within or very close to the baseline region. The 'J35' BD 
sequence and the Skyrmes are also consistent with the MSL baseline results, but the 'PNM' 
BD sequence gives significantly lower pressures than the equivalent MSL sequence. These 
results are expected; the 'J35' MSL and BD EoSs match closely (Fig. 1), but as stated above, 
the BD EoSs give a consistently lower PNM pressure at intermediate densities than the MSL 
EoSs when constrained by the same low density PNM EoS. We can therefore expect any 
quantities that are sensitive to the pressure of the dripped neutrons will also be sensitive to 
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FIG. 17: (Color online) Crust-core and spherical nuclei-pasta transition densities (left) and pres- 
sures (right) versus L for the 'PNM' (thick lines) and 'J35' (thin lines) sequences of the BD EoS, 
together with the selection of Skyrmes. The baseline results are shown for reference as the shaded 
regions. 



these differences in the EoS. 

This is exemplified in Fig. 18, where the pressure scaled by n^^'^ is plotted for the BD EoSs 
on the left and for a selection of Skyrmes on the right. In this figure, and Figs. 19-20, we plot 
both the L = 25, 70 MeV members of the J35 and PNM sequences, the same members that 
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FIG. 18: (Color online) Pressure of matter scaled by n^/^ versus baryon number density in the 
inner crust. On the left the results of the BD EoS are displayed for the L = 25, 35 and 70, 90 MeV 
members of the 'J35' and 'PNM' sequences. The right panel displays the results for our selection 
of Skyrmes together with the DH EoS [38] (DH) and BBP [29] EoSs. The baseline region is shown 
for reference. 
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FIG. 19: (Color online) The volume fraction v occupied by nuclei/nuclear clusters (top), and 
density fraction occupied by free neutrons (bottom), versus baryon number density in the 
crust. The EoSs displayed in the left and right hand panels are the same as Fig. 18. 



largely determine the shaded bounds from the MSL baseline models, and also the L = 35, 90 
MeV members of the same sequences; these bound the set of BD EoSs that are consistent 
with the PNM constraints and 25 < J < 35 MeV. The BD EoSs match the MSL EoSs for 
the low density part of the inner crust, up to n ~ 0.02 fm~^ before diverging significantly 
in the high density region, with the pressure from the BD EoSs being lower, especially for 
lower values of L. However, the Skyrmes used all fall within the baseline region. We also 
plot here, and in Figs. 19-20, the results from the widely used Douchin and Haensel [38] 
(DH) and Baym, Bethe, Pethick [29] (BBP) EoSs, both of which have pressures within the 
baseline region. 

In Fig. 19, V and are plotted versus density for the BD EoSs (left) and the Skyrme 
EoSs together with DH and BBP (right). These quantities are also sensitive to the pressure 
of PNM over the whole density range up to saturation; a lower PNM pressure gives rise to 
nuclei occupying a larger fraction of the space (higher v) and thus a lower density fraction of 
dripped neutrons, which is what we see for most of the crust. Again, the Skyrmes tested are 
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FIG. 20: (Color online) Local proton fraction of clusters x (top) and WS cell size rc (bottom) 
versus baryon number density in the crust. The EoSs displayed in the left and right hand panels 
are the same as Fig. 18. 



consistent with the baseline results. The DH EoS is also consistent, but the BBP EoS gives a 
much lower volume fraction and consequently much higher dripped neutron density fraction 
over the higher density region of the crust. This results from the high compressibility Kq 
used for the BBP EoS and an unrealistically high surface tension. 

In Fig. 20, X and rc are plotted versus density for the BD EoSs (left) and the Skyrme EoSs 
together with DH and BBP (right). The WS cell size rc appears insensitive to the reduced 
neutron pressure of the BD EoSs, while the local proton fraction x drops by up to 0.05 for 
L = 25, 35 MeV; in order to reduce the lattice Coulomb energy cost of increasing the nuclear 
size (which happens because of the reduced dripped neutron pressure), the proton fraction 
drops, which will maintain a roughly constant nuclear charge number. The Skyrmes again 
fall within the region of the MSL baseline results. The DH EoS is also consistent although 
it falls at the upper boundary of the WS cell size region. The BBP EoS predicts a larger 
WS cell size. 
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FIG. 21: (Color online) Crust-core densities scaled by 30 MeV/J (left) versus L and transition 
pressures versus L{n = 0.1 fm~^) - 0.343 i^sym('^ = O.lfm"'^) for the 'PNM' and selected 'J' 
sequences of MSL and BD EoSs, and for selected Skyrmes. 

IX. DISCUSSION 



The transition densities for the sequence of constant J EoSs (Fig. 5) indicates a roughly 
inverse linear scaling of the densities with J. With this in mind, we plot in the left panel of 
Fig. 21 the transition densities multiplied by 30/ J (using J = 30 MeV as a reference value). 
To add the 'PNM' results we must of course use the correlation of J with L. Our Skyrme 
selection is also plotted. We also add onto the plot the BD EoS results; for the BD 'PNM' 
sequence we must also scale the density according to the offset of the correlation relative to 
the MSL correlation. Once this is done, the results form a relatively tight band correlated 
with L. 

It was suggested that a tighter correlation would be obtained by plotting Pec against some 
combination of EoS parameters rather than just L an optimum fit was obtained for 
L{n = O.lfm"^) - 0.343 Ksym{n = O.lfm"'^) using transition pressures obtained using stability 
analyses. Such a correlation can be understood as follows: the pressure of the matter at the 
crust-core transition is related to the pressure of PNM at around the same density (~ 0.1 



fm ); equation 14 shows that this will be related, to the first two orders, with L and 



sym 



at n ~ 0.1 fm~'^ when Kq is fixed. We plot the same quantities obtained in our model for 
the 'PNM' and constant J sequences of both MSL and BD EoSs, as well as the Skyrme 
EoSs, on the right panel of Fig. 21. A tight correlation is indeed observed, reinforcing the 
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FIG. 22: Fraction of mass contained within the pasta phases as a whole and the bubble phases 
relative to the total crust mass versus L. The top left panel shows the baseline results; top right the 
comparison with the BD EoSs and Skyrmes; bottom left and right the variation with the surface 
parameter p. 



applicability of this method, and indicating that more experimental data on the symmetry 
energy at sub-saturation densities, as well as experimental probes of the curvature of the 
symmetry energy i^sym would improve our estimate of the transition pressure. 

The regions over which the MSL and BD EoSs are consistent with both the low density 
PNM constraint and 25 < J < 35 MeV are depicted on the plots; these are still relatively 
large allowed regions. For MSL, the region spans transition densities ~ 0.076 — 0.12 fm~^ 
and transition pressures ~ 0.36 — 0.8 MeV fm~^; for BD the regions are ~ 0.07 — 0.11 fm~^ 
and ~ 0.2 — 0.75 MeV fm~^ respectively. The Skyrmes all fall within the MSL range. 

The transition pressures can be used to estimate the mass and moment of inertia fractions 
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of various components of the crust relative to the crust totals. Assuming an outer layer of 
the neutron star to be thin and contain relatively little mass compared to the whole star, 
the mass of that layer is proportional to the pressure at the base of that layer AM oc P 
|74j . otherwise depending only on the bulk properties of the star. From this, the fraction of 
mass contained within the pasta phases, and within only the bubble pasta phases, can be 
estimated as 

'^^pasta . _ Pp^ A ^bubble . _ Pbu_ ,r,o\ 

7\M ~ P~' 7\M ~ ~P~' ^ ^ 

crust ^ cc crust cc 

where Pbu is the pressure of transition to the bubble phases of pasta. Since the moment of 
inertia contained within the various phases relative to the total crustal moment of inertia is 
proportional to AM, the moment of inertia fractions will be identical to the mass fractions. 

We plot the mass fractions in Fig. 22. In the top left panel, the baseline MSL results are 
shown. For L < 70 MeV, the pasta phases account for a large fraction of the mass of the 
crust: between 50% and 70%. Although the pasta layers occupy a relatively thin layer in 
terms of thickness, the density of those layers is the highest in the crust. For L > 70 MeV, 
the mass fraction levels out for the 'PNM' sequence, always above 40%, while for the 'J35' 
sequence it falls off rapidly to just a few percent at L = 115 MeV. 

Even the bubble phases of pasta can occupy a substantial portion of the crust; between 
10% and 25% at low L < 70 MeV, remaining above 10 % for the 'PNM' seqquence and 
dropping to zero at L 85 MeV for 'J35'. This is of relevance, as in the bubble phases the 
protons are no longer localized in space, opening up a range of possible effects such as the 
direct Urea process ^ and bulk superflow of protons. 

In the top right panel, the MSL baseline is compared to the BD EoS and the selected 
Skyrmes. The mass fractions of the total pasta and bubble phases are elevated relative to 
the MSL results; the total pasta fraction remaining above 60% for the BD J35 sequence, 
and above 40% for the 'PNM' sequence, while the BD 'J35' sequence gives rise to a crust in 
which more than 20% of the mass is contained within the bubble phases. It is important to 
note that this is consistent with the fact that the transition pressures are much lower for the 
BD EoSs; the transition pressures set the absolute mass of the crust, which is lower for the 
BD models; however, the fraction of that mass in pasta and bubbles is higher. The Skyrmes 
have mass fractions in between the MSL and BD EoSs. 

It was noted that the transition pressures are most affected by variations in the behavior of 
the surface energy at low proton fraction, as explored through the variation of the parameter 
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p. In the bottom left and right panels, the baseline results are compared with the results 
using p = 2 (stiffer surface energy at low proton fraction) and p = 4 (softer). For the p=4, 
the mass fractions of pasta and the bubble components are elevated slightly compared to 
the baseline. For p=2 the mass fractions are greatly reduced; remaining between 10% and 
30% for the total pasta fraction on the 'PNM' sequence, and falling from 60% at L =25 MeV 
down to 0% at L =80 MeV on the 'J35' sequence. The bubble phases almost completely 
disappear for p = 2, existing only for L < 50 MeV for any sequence; the bubble fractions 
are raised for p = 4, accounting for 10%-35% of the crustal mass. 

We now give some simple examples of how some bulk crust quantities of interest in 
physical models of neutron star phenomena are affected by the previous exploration of the 
CLDM. We will focus on the charge and mass number of nuclei, the shear modulus of the 
crust and its melting temperature, displayed in Fig. 23. All these quantities are defined only 
within the region of the crust in which separate nuclei exist and become undefined or the 
expression we use is invalid in the pasta phases. 

The top left and right panels show Z and A over the inner crust. The baseline results 
are the shaded band; for a given EoS Z tends to remain roughly constant, and varies by a 
factor of 2 from 20 to 40. Increasing the surface tension at high and low proton fractions 
tends to increase the nuclear size to reduce the surface energy compared to the bulk; thus 
Z increases. The hatched band takes into account those variations through c and p and 
increases Z up to about 50 in the higher density regions of the crust. The DH EoS lies along 
the upper bound of the baseline results; BBP diverges to much higher Z at higher densities, 
again an artifact of the much higher surface tension. A shows a similar trend; the baseline 
results predict A to rise from between 60 and 150 at the lowest densities; higher L predicts a 
slower rise. Increasing the surface energy increases A and the rate at which it increases up, 
reaching towards 1000 at the transition to pasta. Again, DH runs along the upper bound of 
the baseline results, while BBP predicts much higher A at high densities, in line with our 
results for high surface energy. 

The shear modulus of a Coulomb lattice of positively charged nuclei in a uniform neg- 
atively charged background in the neutron star crust at a baryon number density rib was 
determined through Monte-Carlo simulation [2T1423] and can be written as 

/X = 0.1106 (^Y\-'/X^\l-X^f'{Zey, (24) 
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FIG. 23: (Color online) Top left and top right: charge Z and mass A number of nuclei in the 
crust up to the transition to pasta versus density for the baseline MSL EoSs (shaded area) and 
the variation with surface energy (hatched area) together with the results from the Douchin and 
Haensel EoS \^\\ (DH) and BBP EoS [29j (BBP). Bottom left and right: shear modulus of crust 
scaled by pressure and melting temperature versus density for the baseline MSL EoSs. 

where the nuclei are characterized by the nucleon and proton number A, Z and is the 
fraction of neutrons not confined to the nuclei. This is the zero-temperature expression, a 
good approximation for neutron star temperatures below ~ lO^K. Technically this is valid 
only when the ions (nuclei) can be treated as point charges in a uniform background of elec- 
trons; finite size effects and the dripped neutrons are not taken into account. Nevertheless, 
it is widely used to approximate the shear modulus of the inner crust. 

The melting temperature (at which the crystalline lattice changes to a gas of ions) can 
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be expressed as [82] 



neglecting quantum effects lik;e zero-point vibrations. 

The bottom two panels of Fig. 23 show the shear modulus scaled to pressure (left) and 
the melting temperature (right) in terms of 0.01 MeV (1.16 x 10*^ K). The baseline results 
for the shear modulus show a variation of a factor of about 2 throughout the inner crust; 
the melting temperature varies by a factor of 3 throughout the crust. The shear modulus is 
not as sensitive to the surface energy. 



X. CONCLUSIONS 



We have constructed a set of baseline inner crust EoSs using the CLDM, spanning the 
experimentally constrained range of the slope of the symmetry energy 25 < L < 115 MeV 
at constant J over the range 25 - 35 MeV, as well as a sequence in which J is allowed to 
vary in order for the low density EoSs to match current state-of-the-art low density PNM 
calculations. Consistency of the PNM constraint with the range J = 25-35 MeV constrains 
L <70 MeV; the sequences of EoSs consistent with these conditions give the following ranges 
for crustal parameters: 

• rice = 0.08-0.12 fm^^; Up = 0.05-0.06 fm^^. Pasta phases exist to some extent through- 
out the whole range 25 < L < 115 MeV. By scaling rice to the symmetry energy J, the 
correlations with L fall closely along the same curve, independent of the EoS used. 

• Pec = 0.35-0.8 MeV fm'^; Pp = 0.15-0.3 MeV fm'^ We have extended the results 
of [13] to show within the CLDM that P^c correlates tightly with L{n = 0.1 fm~^)- 
0.343i^sym(^ = 0.1 fm~^) independent of J and EoS used. 

• The pressure throughout the crust varies by 5 < P/n^/^ < 15 MeV fm. 

• The volume fraction of the charged component of the crust v ~ 0.001 - 0.15 in the 
region n = 0.001 up to the transition to the pasta phases Up, then rising to the crust- 
core transition volume fraction of 0.8-0.96. 

• Xn rises rapidly from at the neutron drip point to between 0.7-0.9 in the range 
n = 0.005 — Up fm~^; it then falls to its transition value of between 0.01 and 0.2. 

• rc can be roughly be described by the relation rc = —14.5 logri — 3.5±4 fm throughout 
the crust. 
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• X can roughly be described by the relation x — —2n + 0.24 ± 0.03 fm in the range 
n — 0.01 — 0.1 fm"^; at lower densities x rises to 0.3 ± 0.05 at neutron drip density. 

• The fraction of the total crust mass containing the total pasta phases he in the range 
50 - 70%. The mass fractions of the bubble phases, in which the protons are free, is 
in the range 10 - 25%. 

• The charge and mass numbers Z and A fall in the ranges Z = 20-35 throughout most 
of the crust, with the lower bound dropping to Z=10 at the transition to pasta, and 
A = 60-200 throughout most of the crust, with the higher bound rising to 300 at the 
transition to pasta. 

The above ranges are greater when one relaxes the requirement of L < 70 MeV. Note that 
whilst the above ranges give a good guide to the range of crustal parameters, the dependence 
on the symmetry energy within those ranges is more complicated, and to probe the effect 
of the symmetry energy and PNM EoS on neutron star phenomena, one must use crustal 
EoSs consistent with the core EoSs based on the symmetry energy parameters. 

Other model parameters were varied to assess their effect. The surface energy at low 
isospin asymmetry {S — > 0) is increased by decreasing the surface symmetry energy as 
through the correlation with the bulk symmetry energy c. At high asymmetry, the parameter 
p controls the strength of the surface tension, with lower p resulting in a higher surface 
tension. Recent calculations of neutron drops suggest that p might be lower than is obtained 
from fits to calculations involving Skyrme functionals. We show that decreasing p from 3 
(the baseline) to 2 significantly reduces the basehne rice range to 0.06 - 0.11 fm~^ and Pec 
to 0.2 - 0.75 MeV fm~^. The mass fraction of pasta in the crust is reduced to 0.1 - 0.6 for 
L < 70 MeV, vanishing completely for L > 70 in the models with constant J. The effect 
of the surface symmetry-bulk symmetry parameter c is much less pronounced; decreasing c 
from the baseline value to the minimum value used, 2, reduces rice by ~ 0.005 fm^"^ and Pec 
by ^ 0.05 MeV fm~^. Increasing (decreasing) the surface tension can alter Z by up to -|-10 
(-5) throughout the crust. 

The composition is relatively unaffected by changes to c and p except for the WS cell 
size rc which increases by ~ 5 fm by decreasing either c or p to their minimum values. The 
saturation density uq has little effect on any of the crustal properties when the densities are 
scaled to no, except for rc which increases by ~ 5 fm by lowering no to 0.14 fm~^. Changing 
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Kq in the range 220-260 MeV has a small effect on the transition properties, raising or 
lowering n^c by ±0.003 fm~^, with negligible effect on the crustal composition. 

Our results are consistent with other widely used crustal EoSs such as DH. Comparing 
their EoS with the results we obtain using the SLy4 Skyrme (as used in DH), with the 
baseline surface parameters, shows that DH obtain a larger rc (still within our overall 
baseline region) and Z, suggesting that DH used lower values for the surface symmetry 
energy. 

Thus, in order to constrain the crust composition and transition properties further, we 
believe the most urgent work must be done in constraining the PNM EoS, through micro- 
scopic calculations such as those discussed in section HI (GCR, HS, SP) and the surface 
tension of nuclear clusters at low proton fraction. 

We have shown that properties of the crust important in calculating the hydrodynamics 
and thermodynamics of neutron stars (e.g. shear modulus), as well as the thickness, mass 
and moment of inertia fractions of the crust, are sensitive to the EoS of uniform, neutron 
rich matter, especially the symmetry energy (its absolute value at saturation and the shape 
of the symmetry energy function with density), and, relatedly, the PNM EoS, and also to the 
surface tension at low proton fractions. Coupled with the sensitivity of global neutron star 
properties to the same EoS properties, models of phenomena with potentially observable 
signatures such as oscillations, heating and cooling, crustal deformations and the dynamical 
couplings between crust and core can be quite sensitive to the symmetry energy and PNM 
EoS, sensitivities that only appear significant when crust and core descriptions are consistent 
[2]. To that end, the sequences of crust EoSs discussed in this paper are available for use 
[83] in the form of both the basic EoS (pressure versus energy density and baryon number 
density) and an extensive set of crustal composition; appropriate crust EoSs can easily be 
matched to any core EoS based on their respective values of J and L. 

Finally, we have established how the basic crustal composition varies within experimental 
and theoretical uncertainties in the CLDM. The model, as mentioned, has some serious draw 
backs. The WS approximation breaks down in the highest density regions of the crust - the 
pasta layers - and shell effects are not taken into account. In upcoming work we shall 
compare the CLDM results for transition properties and composition with fully microscopic 
calculations of the crustal matter at the highest densities using a 3D Hartree-Fock method 
[54] , in order to establish how much the predictions of the CLDM differ from a more realistic 
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model. 
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XI. APPENDIX: UNIFORM NUCLEAR MATTER EQUATIONS OF STATE 
A. MSL 

The phenomenological modified Skyrme-like (MSL) model |13] will be the form for the 
uniform, isospin-asymmetric nuclear matter EoS we use for our baseline results. The energy 
as a function of density n and isospin asymmetry 5 is written down in a form that closely 
resembles the uniform nuclear matter Skyrme EoS under the Hartree-Fock approach; the 
advantage of the MSL function is that its free parameters can be more easily related to 
the properties of nuclear matter at saturation, and allow for a smooth variation of those 
parameters. The MSL EoS is 



(26) 



^ ' ' n \2ml 2m* ^ ) 

n f n \ ^^^"^ 

E'::Jn) = (1 - y)E%^An,)- + yE':^Jn,) - , (27) 

no V"-o/ 



where 



+n(Ceff + Z^eff5) (28) 



2m* 2m 



2m* 2m 



+ n(Ceff-/^cff5) (29) 



ri = (3/5)(37r^)^/^, uq is the saturation density and a, /3, a, Cqs, D^q, 7sym, -E'sym and y are 
free parameters. Ces and D^h are fixed by setting the effective masses at saturation to be 
m* g = 0.8m and m* q = 0.7m, where m is the average nucleon mass in free space. We set 
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7sym = 4/3; a, /5, a are set by Kq^ uq and Eq while -E'gym and y control the absolute value of 
the symmetry energy and its slope respectively. 

B. BD 

We use as a comparison with MSL a form for the uniform, isospin-asymmetric nuclear 
matter EoS originally obtained by Bludman and Dover [H] and modified by Oyamatsu and 
lida [45]: 

^BD/ S) =-( —rt'l^ + —rt'l^ 



(30) 

n n 

^3 = axn + — , (31) 

1 + 

t^n = h^r? + T^^^. (32) 
1 + O'in 

Vs and are the potential energy densities for SNM and PNM respectively. The free param- 
eters 01,02,03,61,62 are adjusted to obtain the desired values of the saturation properties 
Kq, riQ, Eq, J and L, while 63 which controls the behavior of PNM at high densities is fixed 
at 1.58632 fm^ as suggested in 01 [10]; this parameter is analogous to 7sym in the MSL 
model. 
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